{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "using XSim,JWAS, DataFrames"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "srand(314);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "using XSim"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "100"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "nind     = 100   \n",
    "chrLength= 1.0  #length of each chromosome \n",
    "numChr   = 10   #number of chromosomes\n",
    "nmarkers = 2000 #number of loci for each chromosome\n",
    "nQTL     = 100  #number of QTL for each chromosome"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "build_genome(numChr,chrLength,nmarkers,nQTL)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Sampling 100 animals into base population.\n",
      "Sampling 100 animals into base population.\n"
     ]
    }
   ],
   "source": [
    "popSizeFounder = nind\n",
    "sires = sampleFounders(popSizeFounder);\n",
    "dams  = sampleFounders(popSizeFounder);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Generation     2: sampling    50 males and    50 females\n",
      "Generation     3: sampling    50 males and    50 females\n",
      "Generation     4: sampling    50 males and    50 females\n",
      "Generation     5: sampling    50 males and    50 females\n",
      "Generation     6: sampling    50 males and    50 females\n",
      "Generation     7: sampling    50 males and    50 females\n",
      "Generation     8: sampling    50 males and    50 females\n",
      "Generation     9: sampling    50 males and    50 females\n",
      "Generation    10: sampling    50 males and    50 females\n",
      "Generation    11: sampling    50 males and    50 females\n"
     ]
    }
   ],
   "source": [
    "ngen,popSize = 10,nind\n",
    "sires1,dams1,gen1 = sampleRan(popSize, ngen, sires, dams);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "animals=concatCohorts(sires1,dams1);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "M = getOurGenotypes(animals);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "P = getOurPhenVals(animals,1.0); #residual variance is 1.0\n",
    "nothing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### writeout"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "phenotypes = DataFrame()\n",
    "phenotypes[:y]=P;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "writetable(\"phenotypes.csv\",phenotypes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "writedlm(\"genotype.csv\",M)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## JWAS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "M=readdlm(\"genotype.csv\");\n",
    "phenotypes=readtable(\"phenotypes.csv\");\n",
    "nothing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "20000 markers on 100 individuals were added.\n",
      "MCMC Information:\n",
      "\n",
      "methods                                      BayesC\n",
      "chain_length                                    100\n",
      "starting_value                                false\n",
      "printout_frequency                               20\n",
      "output_samples_frequency                         10\n",
      "constraint                                    false\n",
      "missing_phenotypes                            false\n",
      "update_priors_frequency                           0\n",
      "\n",
      "Information for hyper-parameter: π (Π)\n",
      "π                                              0.95\n",
      "estimatePi                                     true\n",
      "\n",
      "Degree of freedom for hyper-parameters:\n",
      "residual variances:                           4.000\n",
      "iid random effect variances:                  4.000\n",
      "polygenic effect variances:                   4.000\n",
      "marker effect variances:                      4.000\n",
      "\n",
      "\n",
      "\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "running MCMC for BayesC... 14%|████                     |  ETA: 0:00:05"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Posterior means at iteration: 20\n",
      "Residual variance: "
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\r",
      "running MCMC for BayesC... 18%|████                     |  ETA: 0:00:04"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.34606\n",
      "Marker effects variance: 0.271262\n",
      "π: 0.92\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "running MCMC for BayesC... 31%|████████                 |  ETA: 0:00:03"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Posterior means at iteration: 40\n",
      "Residual variance: 0.338365\n",
      "Marker effects variance: 0.35313\n",
      "π: 0.929\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "running MCMC for BayesC... 55%|██████████████           |  ETA: 0:00:02"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Posterior means at iteration: 60\n",
      "Residual variance: 0.34726\n",
      "Marker effects variance: 0.410937\n",
      "π: 0.934\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "running MCMC for BayesC... 75%|███████████████████      |  ETA: 0:00:01"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Posterior means at iteration: 80\n",
      "Residual variance: 0.378563\n",
      "Marker effects variance: 0.47277\n",
      "π: 0.938\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "running MCMC for BayesC... 95%|████████████████████████ |  ETA: 0:00:00"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Posterior means at iteration: 100\n",
      "Residual variance: 0.50035\n",
      "Marker effects variance: 0.53842\n",
      "π: 0.941\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\r",
      "running MCMC for BayesC... 99%|█████████████████████████|  ETA: 0:00:00\r",
      "running MCMC for BayesC...100%|█████████████████████████| Time: 0:00:03\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " 10.836132 seconds (13.00 M allocations: 611.979 MB, 2.98% gc time)\n"
     ]
    }
   ],
   "source": [
    "R=1.0\n",
    "model_equations = \"y = intercept\";\n",
    "model = build_model(model_equations,R);\n",
    "\n",
    "G=0.01\n",
    "add_markers(model,M,G,header=false,G_is_marker_variance=true);\n",
    "\n",
    "@time out=runMCMC(model,phenotypes,Pi=0.95,estimatePi=true,chain_length=100,\n",
    "printout_frequency=20,printout_MCMCinfo=true,methods=\"BayesC\",\n",
    "output_samples_frequency=10);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "using JWAS:misc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "A=GWAS(\"MCMC_samples_for_marker_effects.txt\",model,header=false,\n",
    "window_size=10,threshold=0.01);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "using Plots"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Plots.PyPlotBackend()"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pyplot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<img src=\"\" />"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "plot(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Julia 0.5.0",
   "language": "julia",
   "name": "julia-0.5"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "0.5.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
